Load data and libraries
##################
# LOAD LIBRARIES #
##################
library(tidyverse)
library(Seurat)
library(SeuratObject)
library(tidyseurat)
library(hdWGCNA)
library(enrichR)
library(png)
library(cowplot)
library(patchwork)
library(openxlsx)
source("../../bin/spatial_visualization.R")
source("../../bin/plotting_functions.R")
#########
# PATHS #
#########
input_dir <- "../../results/06_DGE_condition_st_data/"
result_dir <- "./Figures/05/"
if( isFALSE(dir.exists(result_dir)) ) { dir.create(result_dir,recursive = TRUE) }
epi_clus <- "^5$|^6$|^7|^8" # res 0.7
ord <- c("Superficial", "Upper IM", "Lower IM", "Basal","1","4","0","3","2","9","10","11","12")
ord1 <- c("5", "6", "7", "8","1","4","0","3","2","9","10","11","12")
sample_id <- c("P020", "P045", "P050", "P057",
"P008", "P031", "P044","P080", "P026", "P105",
"P001", "P004", "P014", "P018", "P087", "P118",
"P021", "P024", "P067", "P081", "P117" ) %>% set_names()
#############
# LOAD DATA #
#############
# hdWGCNA
DATA <- readRDS(paste0("../../results/09_hdWGCNA/","hdWGCNA_3771DEGs_Seurat.RDS"))
modules <- read_csv(paste0("../../results/09_hdWGCNA/", "wgcna_3771DEGs_modules.csv"))
# DATA <- readRDS(paste0("../../results/09_hdWGCNA/",,"all_Clus_4000DEGs/","hdWGCNA_Seurat.RDS"))
# modules <- read_csv(paste0("../../results/09_hdWGCNA/all_Clus_4000DEGs/", "wgcna_all_Clus_modules.csv"))
# plot the dendrogram
png(paste0("./Figures/05/", "Spatial_hdWGCNA_dendrogram.png"),
width = 10, height = 3, units = "in", res = 300)
PlotDendrogram(DATA, main=NA, marAll = c(1, 4, 1, 0))
dev.off()
## quartz_off_screen
## 2
dend <- readPNG("./Figures/05/Spatial_hdWGCNA_dendrogram.png")
g <- grid::rasterGrob(dend, interpolate=TRUE)
A <- ggplot() +
annotation_custom(g, xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=Inf) +
theme_nothing() +
theme(rect = element_blank(), # removes the box around the plot)
plot.margin = unit(c(-0,0,0,0), "lines")) #t,r,b,l
# dev.new(width=3.5, height=3, noRStudioGD = TRUE)
A

# get module eigengenes and gene-module assignment tables
MEs <- DATA@misc[["vis"]][["MEs"]]
# add the MEs to the seurat metadata so we can plot it with Seurat functions
DATA@meta.data <- cbind(DATA@meta.data, MEs)
# plot with Seurat's DotPlot function
mods <- unique(modules$module)[unique(modules$module) != 'grey']
p <- c('groups', 'layers') %>%
map(., ~DotPlot(DATA, features=rev(mods), group.by = .x, dot.min=0.1) )
# flip the x/y axes, rotate the axis labels, and change color scheme:
p <- map(p,~.x +
coord_flip() +
RotatedAxis() +
scale_color_gradientn(colours = c('blue','grey95', 'red'), name = "Avg. Expression",
limits = c(-1.5, 2.5), oob = scales::squish,
values = scales::rescale(c(-1.5, 0, 2.5))) +
theme(axis.text = element_text(size = 10),
legend.text = element_text(size = 9), legend.title =element_text(size = 10),
axis.line = element_line(size = .4), # , colour = "#bebebe"
axis.ticks = element_line(size = .4), # , colour = "#bebebe"
legend.direction = "vertical", legend.box = "horizontal",
legend.margin=margin(1,-1,1,6),
plot.margin = unit(c(.7, -0, -1.2, -1), "lines")) + # t,r,b,l
xlab('') + ylab('') + guides(col = guide_colourbar(barwidth = .3, barheight = 4 ))
)
# combine plots
# dev.new(width=3.9, height=3.3, noRStudioGD = TRUE)
leg <- get_legend(p[[1]])
p1 <- plot_grid(p[[1]]+theme(legend.position = "none"),NULL, rel_widths = c(1, 1.3))
(B <- plot_grid(p1, p[[2]]+theme(legend.position = "none"),rel_heights = c(.75,1,1), ncol = 1) )

ggsave(filename=paste0("./Figures/05/", "dot-plot.png"),B, width = 3.9, height = 3.3, bg = "white")
# relative expression level of each module
plot_list <- ModuleRadarPlot(
arrange(DATA, layers),
group.by = 'groups', combine = F, # ncol = 4,
#barcodes = seurat_obj@meta.data %>% subset(cell_type == 'INH') %>% rownames(),
axis.label.size=3, group.line.width = 0.3, grid.line.width = 0.3, gridline.max.linetype = "dashed",
grid.label.size=3
)
# dev.new(width=5.5, height=2, noRStudioGD = TRUE)
(C <- wrap_plots(plot_list, ncol=4) & theme(title = element_text(size=8), plot.margin = unit(c(.1, -0, -1, -0), "lines")) )

# ggsave(filename=paste0("./Figures/05/", "Group_contribution.png"),C, width = 10, height = 3, bg = "white")
plot_filt.fun <- function(DATA, gr = "L1"){
DAT <- filter(DATA, groups == gr)
DAT@misc[["vis"]][["MEs"]] <- DATA@misc[["vis"]][["MEs"]][colnames(DAT),]
plot_list <- ModuleFeaturePlot(DAT, reduction = "umapharmony", features = "MEs", title =F)
return(plot_list)
}
mod <- c("SM1", "SM2", "SM3", "SM4") # , "SM5", "SM6"
mod <- map(mod, ~plot_genes.fun(DATA,
gene = .x,
scale = F,
mins = -20, maxs = 20,
diverging = T,
col = rev(c("#c41625","#dc4e43","#fa9975","#FDDBC7","#F7F7F7","#D1E5F0","#92C5DE","#4393C3","#2166AC")),
point_size = .5,
red="umapharmony",
lable = TRUE) + ggtitle(" ") + theme(plot.title = element_text(hjust = 0))) # size = 3
# dev.new(width=4, height=4, noRStudioGD = TRUE)
(D_1 <- map(seq_along(mod), ~mod[[.x]] + facet_wrap(~groups, ncol = 4)) %>% wrap_plots(., ncol=1) )

# ggsave(filename=paste0("./Figures/05/", "Modules_across_groups_UMAP.png"),D_1, width = 10, height = 10, bg = "white")
# enrichr databases to test
dbs <- c('GO_Biological_Process_2021','KEGG_2021_Human','Transcription_Factor_PPIs')
# perform enrichment tests
DATA <- RunEnrichr( #map(dbs, ~
DATA,
dbs=dbs, # character vector of enrichr databases to test
max_genes = Inf # number of genes per module to test. use max_genes = Inf to choose all genes!
)
## Uploading data to Enrichr... Done.
## Querying GO_Biological_Process_2021... Done.
## Querying KEGG_2021_Human... Done.
## Querying Transcription_Factor_PPIs... Done.
## Parsing results... Done.
## Uploading data to Enrichr... Done.
## Querying GO_Biological_Process_2021... Done.
## Querying KEGG_2021_Human... Done.
## Querying Transcription_Factor_PPIs... Done.
## Parsing results... Done.
## Uploading data to Enrichr... Done.
## Querying GO_Biological_Process_2021... Done.
## Querying KEGG_2021_Human... Done.
## Querying Transcription_Factor_PPIs... Done.
## Parsing results... Done.
## Uploading data to Enrichr... Done.
## Querying GO_Biological_Process_2021... Done.
## Querying KEGG_2021_Human... Done.
## Querying Transcription_Factor_PPIs... Done.
## Parsing results... Done.
# retrieve the output table
enrich_df <- GetEnrichrTable(DATA) %>%
filter(Adjusted.P.value < 0.05) %>%
split(~db)
enrich_df %>%
map(., ~split(.x, ~module)) %>%
imap(., ~write.xlsx(.x, paste0(result_dir,"New_3771DEGs/", "Enrichment_",.y,".xlsx")) )
## $GO_Biological_Process_2021
## A Workbook object.
##
## Worksheets:
## Sheet 1: "SM1"
##
##
## Sheet 2: "SM2"
##
##
## Sheet 3: "SM3"
##
##
## Sheet 4: "SM4"
##
##
##
## Worksheet write order: 1, 2, 3, 4
## Active Sheet 1: "SM1"
## Position: 1
##
##
## $KEGG_2021_Human
## A Workbook object.
##
## Worksheets:
## Sheet 1: "SM1"
##
##
## Sheet 2: "SM2"
##
##
## Sheet 3: "SM3"
##
##
##
## Worksheet write order: 1, 2, 3
## Active Sheet 1: "SM1"
## Position: 1
##
##
## $Transcription_Factor_PPIs
## A Workbook object.
##
## Worksheets:
## Sheet 1: "SM1"
##
##
## Sheet 2: "SM2"
##
##
## Sheet 3: "SM3"
##
##
## Sheet 4: "SM4"
##
##
##
## Worksheet write order: 1, 2, 3, 4
## Active Sheet 1: "SM1"
## Position: 1
# saveRDS(enrich_df, paste0("../../results/09_hdWGCNA/New_3771DEGs/", "Enrichment.RDS"))
# enrich_df <- readRDS(paste0("../../results/09_hdWGCNA/New_3771DEGs/", "Enrichment.RDS"))
#######################
# ENRICHMENT BARPLOT #
#######################
overlap.fun <- function(string){
l <- str_split(string, pattern ="/")
l <- map_dbl(l, ~as.numeric(.x[1])/as.numeric(.x[2]) )
return(l)}
# str_match("cytoplasmic translation (GO:0002181)", "^(.+?)\\s\\((.+)\\)$")[2]
GeneRatio_plot.fun <- function(enrich_df, txt_size = 15, nr_path=3,
col=c("#ed968c","#f9d14a","#88a0dc","#e78429")){
dot_df <- enrich_df %>% # dot_df <- enrich_df$GO_Biological_Process_2021 %>%
{if(grepl("GO",.$Term[[1]])){mutate(., "Term" = str_match(.$Term, "^(.+?)\\s\\((.+)\\)$")[,2],
"GOid" = str_match(.$Term, "^(.+?)\\s\\((.+)\\)$")[,3] )}else .} %>%
filter(Adjusted.P.value < 0.05) %>%
mutate("-log10(P-value)" = -log10(P.value)) %>%
mutate(GeneRatio = overlap.fun(.$Overlap)*100)
dot_df <- dot_df %>%
group_by(module) %>%
top_n(., -nr_path, Adjusted.P.value)
p <- ggplot(dot_df, aes(x = `-log10(P-value)`, y = fct_reorder(Term, `-log10(P-value)`), fill = module, col = module)) +
geom_col(width = .05, show.legend = F) +
geom_point(aes(size = GeneRatio)) + theme_classic() +
scale_fill_manual(values = col, aesthetics = c("fill", "colour")) +
facet_wrap(~module, scales = 'free', ncol = 1) +
#scale_x_continuous(expand = c(0, 1.5)) +
coord_cartesian(clip = F) +
scale_x_continuous(limits = function(x){c(0, +max(0.1, x))}) +
scale_size_continuous(breaks = c(10, 50, 80)) +
theme(axis.title.y = element_blank(),
strip.text.x = element_text(hjust = 0.1, margin=margin(l=0)),
strip.background = element_blank(),
panel.spacing = unit(1, "lines"),
axis.text.y = element_text(size = 10),
panel.border = element_blank())
return(p)
}
df <- bind_rows(enrich_df[[2]], enrich_df[[3]]) %>%
filter(db == "Transcription_Factor_PPIs" & module == "SM4" | db == "KEGG_2021_Human")
# dev.new(width=3.5, height=5, noRStudioGD = TRUE)
(D_2 <- GeneRatio_plot.fun(df) + theme(legend.position = "none") )

# ggsave(paste0("./Figures/05/", "Enrichment.png"),D_2, width = 3.8, height = 5, bg = "white")
######################
# COMBINE ALL PANELS #
######################
# dev.new(width=7, height=8, noRStudioGD = TRUE)
A_C <- plot_grid(A, C, ncol = 1, rel_widths = c(1, .4))
A_B_C <- plot_grid(A_C, B, ncol = 2)
D <- plot_grid(D_1, D_2)
(Figure5 <- plot_grid(A_B_C, D, ncol = 1, rel_heights = c(0.33, .66)) )
